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ABSTRACT 

We solve the inviscid Euler equations in complicated geometries using a Cartesian 
grid. This requires solid wall boundary conditions in the irregular grid cells near the 
boundary. Since these cells may be orders of magnitude smaller than the regular grid 
cells, stability is a primary concern. We present a.new approach to this problem and 
illustrate its use. 
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1 Introduction 


In previous work [1], [6], [7], we have described a Cartesian grid method for the 
inviscid Euler equations in arbitrary geometries. There are many advantages to be 
gained from this approach. Grid generation is simplified, since we avoid the use of 
(possibly multiblock) body-fitted grids, and we can use high resolution, highly efficient 
solvers on regular grids over the bulk of the domain. This has led to renewed interest 
in Cartesian grids in recent years, e.g., [3], [10]. One of the difficulties with Cartesian 
grids is that they give insufficient resolution in certain regions such as leading edges. 
This can now be overcome by Cartesian adaptive mesh refinement [1], [2], 

The principal remaining difficulty in this approach is due to the essentially ar- 
bitrary way that a Cartesian grid intersects the boundaries of the computational 
domain. In particular, a solid wall boundary cutting through the grid creates irreg- 
ular cells that may be orders of magnitude smaller than the regular cells away from 
the boundary. For these irregular cells, special difference equations are needed that 
maintain stability and accuracy, and satisfy the solid wall boundary conditions of no 
normal flow. 

In this work, we present an improved method for the small boundary cells. We 
use an explicit, finite volume formulation that computes fluxes at cell edges on the 
regular part of the domain. We would like to define fluxes at the edges of the irregular 
cells in such a way that the method is stable even with very small boundary cells, 
using a time step based on the regular grid cells away from the boundary. The CFL 
condition requires that the numerical method allow information to propagate at least 
as quickly as the underlying differential equation. In the present context this means 
that we must define fluxes at the sides of our irregular cells based on more than just 
the neighboring cell values. 

In our previous work, we have used a wave propagation approach in defining these 
fluxes. Here we propose an alternative method that has some advantages over the 
wave propagation approach. In particular, the wave propagation method is subject 
to intermittent instabilities due to two-dimensional effects that are not clearly under- 
stood. The new method has a cancellation property in two dimensions that appears to 
give better stability properties. Moreover, the computational geometry is simplified 
in the new approach. The fluxes are defined in terms of weighted averages of nearby 
cell values. These weights may be calculated as a preprocessing step on any fixed grid 
and need not be repeatedly calculated. In the previous approach the weights depend 
on the flow variables and a certain amount of computational geometry was required 
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near the boundary in every time step. 

We consider the inviscid Euler equations of gas dynamics in two space dimensions, 

+ /(«)* + 9{u)v = 0 ( 1 . 1 ) 

where 
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Here (111,112) represents the velocity, E is the total energy per unit mass, and p is the 
pressure, which is related to the other variables by the equation of state. We assume 
a 7-law gas, so that 

P = {l -WpE -^p{u\+u\)). ( 1 . 3 ) 

At a solid wall boundary we require that the component of velocity normal to the 
wall be zero. 

In one space dimension the system reduces to 

u t + f(u) x = 0 ( 1 . 4 ) 

where u = ( p , pv, pE ) and f(u) = ( pv,pv 2 + p,v(pE + p)), with v = u x the velocity. 
The boundary conditions become v = 0 at a solid wall. 

2 A one- dimensional example 

In order to illustrate this approach we begin with a one-dimensional model problem, 
the one-dimensional Euler equations for x > 0 with a solid wall at x = 0 . We take a 
grid with cell interfaces at the points 

Xq = 0 

X\ = h' 

Xj = h! + jh for j = 2 , 3 , 

Here h is a uniform grid spacing and h! < h. the grid is uniform except for one small 
cell near the boundary (see Figure 1 ). We use a conservative method in the form 

o? +1 = v? - - q 1 ], i = 0, 1, ... . ( 2 . 5 ) 
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Figure 1: One dimensional grid with one irregular cell adjacent to the wail. 

Here hj is the width of the jth cell, so in our case we have ho = h' and hj = h for 

j > 0. 

For simplicity we restrict our attention to Godunov’s method on the regular por- 
tion of the grid, although the ideas we propose can be extended to higher order 
methods as well. In Godunov’s method we take 

*?-/(»•(£$-.,£?)) P- 6 ) 

where u*(u L ,u R ) represents the solution to the Riemann problem with left and right 
states u L and tt*, evaluated along x/t = 0. Although a rigorous stability proof is not 
available for systems of equations, in practice this method is always stable provided 
the CFL condition 

< 1 (2.7) 

is satisfied, where is the maximum wave speed. We will assume that our time 
step k is chosen so that the condition (2.7) is satisfied relative to the uniform h. We 
will use the flux (2.6) for j = 2, 3, . . ., i.e., at all interfaces where the cell on both 
sides is regular. Our task is to define fluxes F J* for j = 0, 1 so that we maintain 
stability (and accuracy) with this time step even if h ' -C h. 

First suppose h' = h. Then we can use the Godunov flux (2.6) also at j = 1. 
At the wall we use the well-known observation that the solution to the boundary 
value problem can be obtained by ignoring the wall and extending the computational 
domain to the whole line — oo < x < oo, if we take data uo(x) for x < 0 equal to 

/>(x,0) = p(-x,0) 

u(x,0) = — u(— x,0) for x < 0 

p(x, 0) = p(— x,0). 

We will denote this “reflection” of the data (in which the velocity is negated) by the 
operator 71, so that for shorthand we can write 

u(x,0) = 7£(u(-x,0)) for x < 0. 
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With this extended data, the solution continues to satisfy u{x,t) = Tl(u[— x,t)) also 
for t > 0 and in particular the boundary condition u(0,t) = 0 is automatically 
satisfied. This suggests that we obtain a flux at the wall by solving a Riemann 
problem with left and right states 

u L = TZ(U 0 ), u R = U 0 

in each time step to obtain 


F 0 = f(u*(Tl(Uo),uo). 


(For brevity we will leave off the superscript n in general.) Note that the density 
and energy components of this flux will be zero since the velocity component of u* is 
zero at the wall. There will only be a momentum flux at the wall due to the pressure 
there, as expected physically. 

If h' <. h we could attempt to use this same formula to define Fq but we would 
find that it is unstable unless the CFL condition 


kK 


h' 


< 1 


( 2 . 8 ) 


is satisfied. This will place an unreasonable restriction on k if h' <C h. 

This instability is caused by the fact that the boundary flux F 0 is based on the data 
U 0 alone. If the CFL condition (2.8) is satisfied, then it is only this data that affects 
the flux at the wall over the time step. However, when (2.8) is violated the value Uy 
should also affect the flux at the wall, and ignoring this effect leads to instability. 

In a “large time step” approach we increase the stencil of the method, meaning we 
allow more data points to come into the computation of each flux, and hence retain 
stability. One way to achieve this is by a wave propagation approach. The solution of 
the Riemann problem at each cell interface consists of three waves propagating away 
from the interface. If (2.8) is satisfied then these waves remain in the cells bordering 
the interface during the entire time step and hence affect the solution only in these 
cells. If (2.8) is violated then the waves may affect cells further away. Implementing 
Godunov’s method in terms of this wave propagation approach, allowing waves to 
affect more than just the adjacent cell, gives a large time step generalization that 
remains stable for much larger time steps[5]. In the present context this allows us 
to reduce h' without reducing the time step k. Waves from the boundary Riemann 
problem cross the interface at Xy and affect Uy as well as U 0 . Waves from the interface 
at Xy may reach the boundary. These waves reflect off the boundary and the reflected 
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wave affects the value U 0 and perhaps also Ui if the reflected wave reaches the cell 
interface at x\ during the time step. 

A more detailed description of this procedure may be found in [7]. A natural 
extension to two space dimensions gives one method to deal with small cells near 
the boundary, as described in [1], [6], [7]. In one dimension this works very well 
but in two dimensions occasional stability problems have still been observed due to 
multidimensional effects. 

The new approach. Our new approach to the small cell problem can also be 
illustrated with the one-dimensional problem described above. We again use the 
method (2.5) with Godunov fluxes (2.6) for j = 2, 3, — For j = 0 and j = 1 we 
define fluxes in a similar manner but with a different choice of states u L and u in the 
Riemann solver. Recall that in a naive attempt to use Godunov’s method regardless 
of the size h' of the small cell we would take left and right states 

u L 0 = K(U 0 ) < = U 0 (2-9) 

u\ = U 0 uf = Ul (2.10) 

To maintain stability when b! is small, we need to allow data from additional grid 
cells to affect the left and right states at each of these interfaces. Recall that the 
method is assumed to be stable with our choice of k and h on the regular portion of 
the grid. This suggests that we should define Uj by taking the average value of U 
over an interval of length h to the left of the interface Xj and define u* by taking the 
average value of U over an interval of length h to the right of Xj. 

For example, at Xo = 0 (the wall) we set 

uZ = Uh'U 0 + (h- h!)lh) (2.11) 

ft 

If we view the grid values as defining a piecewise constant function with values Uj 
in the jth cell, then (2.11) is the average value of this function over the interval 
0 < x < h. Note that if h' = h (the grid is Completely regular) then (2.11) reduces 
to Uq = U 0 as expected for Godunov’s method. Recall that in Godunov’s method 
we take Uq = 7l(Uo) = 7£(uJ) to impose the boundary condition u(0,t) = 0. This 
suggests that more generally we take 

uS = ft«) (2.12) 

where Uq is defined by (2.11). We then use the Godunov flux 

fb = /(«•(«$, O) ( 2 - 13 ) 
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as the flux at the wall. Using (2.12) guarantees that there will be no flux of mass or 
energy through the wall and hence that the method is conservative. 

To define the left and right states at xi we again construct intervals of length h 
to either side of this point and average the piecewise constant function defined by U 
over these intervals. To the right oi x\ lies a regular cell of length h and so 

u? = Ui. (2.14) 

To the left of X\ an interval of length h extends beyond the wall (assuming h' < h). 
Beyond the wall we assume that U takes the value u% given by (2.12). A weighted 
average of this value and Uo gives u\\ 

u\ = \{h'U 0 + {h- fc'K). (2.15) 

The flux fi is then defined by 

= ( 216 ) 

Again, if h' = h this reduces to the standard Godunov flux. 

This method remains stable even when h' <C h. To see why this should be so, 
consider the formula (2.5) for j = 0 where hj = h ' . It is the division by h! that 
may cause stability problems unless the fluxes Fq and F\ themselves agree to 0(h') 
as h' — > 0. The Godunov fluxes based on (2.9), (2.10) do not have this property. 
However, our proposed fluxes (2.13) and (2.16) do have this property, since inspection 
of the formulas (2.11), (2.12), (2.14), and (2.15) shows that u\ = u% + O(h') and 
Uq = U* -f O(h') as h' — ♦ 0. Since the flux function /( u*(u L ,u R )) is a Lipschitz 
continuous function of u L and u R , it follows that Fi — F 0 — O(h') as h' —* 0 and there 
is at least a chance that the method remains stable for arbitrary h' <C h. Numerical 
experiments show that this is indeed the case (although it is possible to contrive 
examples, such as a strong rarefaction wave originating at this irregularity, where the 
results are not very accurate). 

3 Boundary conditions in two dimensions 

Turning now to the two-dimensional problem, we will give a brief description of how 
the idea described above extends to handle the small cell problem. 

Consider the portion of the boundary shown in Figure 2a and a typical boundary 
cell (i,j). The formula for updating the value Uij is the two-dimensional analog of 

(2.5). 

tf? 1 = V it - ~ Fa + - G„ + H it \. (3.17) 
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Figure 2: (a) The Cartesian grid near the boundary, (b) Blow up of cell (i j) showing 
the location of fluxes. 

The fluxes F, G, and H represent flux per unit time through the corresponding side 
of the grid cell (see Figure 2b) and Aij is the area of the cell. If any of the sides are 
missing, the corresponding flux is zero. 

On regular grid cells, H%j = 0 and the fluxes F and G might be defined by an 
extension of the Godunov method, setting 

Gii = hgWij-uUii)). (3.18) 

Here u* represents the solution to the appropriate one- dimensional Riemann problem 
in the x or y direction. Note that the fluxes include the factor h, the length of each 
side, to give a flux per unit time across the side. 

It is the denominator Aij in (3.17) that causes trouble when the cell is very small. 
We again assume the method is stable on the regular portion of the grid, where 
Aij = h 2 . To maintain stability we need to insure that our formulas for the fluxes 
cause the total flux (the sum in brackets in (3.17)) to cancel to O(Ajj) as Ajj — » 0. 
This is only possible if the fluxes axe computed via formulas that involve more than 
just the two cells bordering the cell side. We take an approach analogous to what we 
described above in one dimension. 

Boundary fluxes. We begin by considering the boundary segment, where we 
must compute the flux Hij. In two dimensions the solid wall boundary condition 
requires that the normal velocity at the wall be equal to zero. If we have some value 
u'jj representing the value of u just inside the wall, then we can obtain the flux Hij 
by solving a one- dimensioned Riemann problem in the direction normal to the wall, 
with left and right states 

< = tig- = «§'• 
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Figure 3: The inbox and outbox constructed from the boundary segment of cell ( i , j), z 

and the inbox for two neighboring cells. % 

The reflection operator 71 is now defined by negating the normal velocity component 
while leaving the tangential velocity component along with the density and pressure 
unchanged. The resulting Godunov flux is used for Hij. 

We obtain it- by a procedure analogous to the one-dimensional example. We 
construct a box extending a distance h away from the wall as shown in Figure 3. The 
box extending into the computational domain is called inbox(i, j). The mirror image 
box outside the domain is called outbox(i, j). We obtain the value u £ by viewing 
the given data U as defining a piecewise constant function, constant in each grid cell, 
and setting u £ to be the average value of this function over the region inbox(i, j). In 
Figure 3 inbox(i, j) would contain an area- weighted average of two grid values while 
the value for inbox(t + l,j) is based on four grid values. We think of the outbox as 
containing the value u-' = 

To find the weights needed to compute we must compute the intersection of the 
inbox with each nearby cell. This is easily accomplished with standard computational 
geometry routines. Note that for a given geometry and grid these weights need only 
be computed once at the beginning of the computation. They need not be recomputed 
in each time step. 

Fluxes at other sides. We now consider the fluxes F and G along other sides 
of this cell. These are all computed by similar procedures, so to be specific we will ~ 

consider the computation of Fy, the flux on the left side of this cell. 

To compute Fij we solve two Riemann problems, one in some direction £ with 
some data uf, u* and the other in the orthogonal direction r] with data u%, u*. The - 

choice of these directions and the data will be discussed in a moment. First we explain 
how these Riemann problem solutions are computed and used to define F,j. 

Figure 4 shows a typical vertical cell interface and two orthogonal directions <f and 
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Figure 4: A vertical cell interface and the £- and 77- direct ions. 

rj. Let 9 be the angle that £ is rotated from the x-direction (9 < 0 in this example). 
Suppose we solve a one- dimensional Riemann problem in the ^-direction with left and 
right states u £, u* to obtain the flux per unit length per unit time in the ^-direction. 
(To do this we rotate the velocity components of u|, u* into £-77 velocity components, 
solve the one- dimensioned Riemann problem, and then rotate the resulting flux / back 
to x-y velocity components.) Call this resulting flux /$. 

Similarly, solving a one-dimensional Riemann problem in the 77-direction with 
left and right states u%, u* gives f v , the flux per unit length per unit time in the 
77-direction. The total flux across the vertical segment of length h' is then 

F = h\f( cos 9 — f v sin 9). (3.19) 

This is the value we use for the flux F t j . 

This same approach has been used by others ( e.g ., [4], [8], [9]) to define multi- 
dimensional upwind methods. In these methods the directions ( and 77 are chosen 
based on the local flow in an attempt to use physically meaningful directions in place 
of the artificial coordinate directions. For example, the direction of the velocity or 
the pressure gradient might be used to define (. In our application we are only con- 
sidering cells adjacent to the boundary and the relevant directions are the directions 
tangential and normal to the wall. We choose £ to be the direction tangential to the 
wall in one of the two cells bordering this interface. Since our primary concern is to 
maintain stability in very small cells, we choose the smaller of the two adjacent cells 
to define this direction. This will lead to cancellation of fluxes in tiny cells in the 
same manner as previously seen in the one-dimensional example. The 77-direction is 
normal to the ^-direction. 

Tangential boxes. We must still specify the data for these tangential and normal 
Riemann problems. We first consider the tangential problem. We use an approach 
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Figure 5: (a) Tangential boxes constructed from the cell interface, (b) Normal boxes 
from the cell interface and outbox(t — l,j). 

similar to the specification of data in an inbox described above. From the interface 
we construct boxes that extend a distance h in the ^-direction. Figure 5 a shows an 
example. The data u%, u* is obtained by an area- weighted average of the values 
in each cell the box overlaps. In our current implementation we assume the wall 
is convex, so that these boxes lie entirely within the computational domain. Each 
box overlaps at most two grid cells and the weights are easily calculated. Since the 
directions £ and rj and the resulting boxes depend only on the geometry, not on the 
flow variables, these weights can again be calculated once and for all as a preprocessing 
step. 

Normal boxes. Figure 5 b shows the normal boxes in the 77 -direction. The box 
in the outward direction does not hit the boundary and overlaps at most two regular 
cells, so u* is calculated as an area- weighted average of these cell values. The other box 
may extend beyond the boundary. If so, the portion lying outside the computational 
domain lies in one or more out boxes, the artificial cells created in the process of 
computing the boundary flux Hij described above. Figure 5b shows a simple example 
where the normal box intersects only one cell (z — 1 ,j) and outbox(z — 1, j). More 
generally the normal box might intersect two cells and their outboxes, as happens 
for example when we compute the flux F i+ x, j which involves cells (i,j) and (z, j — 1 ). 
Moreover the two outboxes will in general overlap due to the convexity oTthe region. 
We again use area-weighted averaging over the four cells in question, weighting the 
values Uij, by the areas of intersection and then dividing by the 

sum of all these areas. 

Cancellation. Although we will not present the details here, it can be shown 
that this way of defining fluxes leads to the desi red cancellation of fluxes in v ery small 
cells. The values computed at each of the three sides of a very small triangular 
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cell are nearly the same because of our construction. They differ by only O ( A XJ ) as 
Aij — * 0. The same is true of each of the other values u*, u^, U* and so by Lipschitz 
continuity of the fluxes F, G and H we obtain the required cancellation. Numerical 
results show stability even when A{j is many orders of magnitude less than h 2 . 

Higher order methods The method (3.17) with fluxes (3.18) is only first order 
accurate and is highly dissipative. In our previous work we used the wave propagation 
boundary conditions together with a high resolution method away from the boundary 
and obtained reasonable results (e.g., [1]). The new boundary conditions can also 
be applied in conjunction with a high resolution method and gives similar results. 
Moreover, with our new formulation it appears to be easier to improve the accuracy 
of the boundary conditions, allowing us to obtain higher order accuracy overall. The 
main idea is to introduce slopes in each cell and use piecewise linear approximations 
in place of piecewise constants to define the fluxes. Near the boundary we can easily 
estimate slopes in the tangential direction along the wall by differencing values in the 
inboxes that we have defined above. 

These improvements are still being investigated and will be reported in detail 
elsewhere. Here we will only compare results obtained with the method as we have 
described it and results obtained using the same interior method with the wave prop- 
agation boundary conditions described in earlier papers. 


4 Numerical results 


We show one representative test case, a supersonic shock going around an expansion 
corner. We also show the steady state solution obtained at large times. The exact 
rarefaction wave solution is a simple wave and can be computed following Section 
6.17 of Whitham[ll], for example. 

The geometry we use is the rectangle [0, 1.32] x [0,0.8] with a solid wall at 


V = 


0.3 

0.3(1 -(x-0.1) 2 ) 
0.192 - 0.36(x - 0.7) 


x < 0.1 

0.1 < x < 0.7 

0.7 < x < 1.32. 


The initial conditions consist of a Mach 2.31 shock at x — 0.06 with left and right 


states 

/= 5.1432, i4 = 2.045U, u\ = 0, p* = 9.04545 


and 

p R = 1.4, uf = 0, u% = 0, p* = 1.0. 
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Figure 6: Shock propagation results at t = 0.4 (a) Using the wave propagation 
boundary conditions, (b) Using the new boundary conditions. 

We take h = 0.02 (66 x 40 grid) and a time step k = 0.002. This corresponds to a 
Courant number of roughly 0.37 relative to the regular cells with area h 2 . For the 
crude form of Godunov’s method used here, the stability restriction requires Courant 
number less than 0.5. The smallest cells near the boundary have an area roughly 

10" 3 /i 2 . 

Figure 6 shows numerical results at time t — 0.4, as the shock is rounding the 
corner. Results obtained with the wave propagation boundary conditions are shown 
in Figure 6a, while Figure 6b shows the results obtained with our new approach. 
These results axe very similar. Slight discrepancies can be seen near the wall just 
aro un d the shock. For this problem both sets of boundary conditions worked well. 
We have also performed tests on other problems where the wave propagation method 
shows instabilities and have observed no such difficulties with the new method. 

Figure 7 shows the steady state results obtained after many iterations of the time 
dependent code (no attempt has been made so far to accelerate convergence for steady 
state solutions). We only show the results with our new boundary conditions. The 
wave-propagation boundary conditions give very similar results. We use a coarser 
grid than in the previous example ( h = 0.4) in order to demonstrate that we achieve 
reasonable accuracy along the boundary even with a relatively coarse piecewise lin- 
ear representation of the boundary. We also use a larger computational domain, 
[0,2] x [0, 1.6] to minimize the impact of the far-field boundaries. The true solution 
is a rarefaction wave originating from the portion of the boundary with nonzero cur- 
vature. In the exact solution the contour lines would be straight lines. Our results 
are contaminated by effects from the far-field boundary. 

Near the solid wall the contour lines appear to show a boundary layer. This is an 
artifact of the graphics routine, which assumes the data is on a uniform grid at cell 
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Figure 7: Steady state results, (a) Pressure contours, (b) Pressure along the wall. 
The solid line is the exact solution. + marks are the numerical solution. 

centers. Our data near the boundary should be viewed as an approximation to the 
pointwise value at the center of mass of the irregular cell, not at the center of the full 
Cartesian cell. 

In order to examine the accuracy at the wall, Figure 7b shows plots of the pressure 
along the boundary, plotted against arclength. To obtain a boundary pressure, the 
cell value U {j and the reflected value 7 l{Uij) are used to solve the one-dimensional 
Riemann problem normal to the wall in each irregular cell. The resulting pressure p 
is used as the boundary pressure. Figure 7b shows these results and also the exact 
solution (to machine precision) calculated using the theory of [11]. 

In more complicated computations we use adaptive grid refinement to obtain high 
resolution results with minimal effort. The boundary conditions described here can 
also be used in conjunction with the adaptive Cartesian grid code described in [1] 

and [2]. 
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